Assignment: Notebook for Peer Assignment

Estimated time needed: 60 minutes

Assignment Scenario

Congratulations! You have just been hired by a US Weather forecast firm as a data scientist.

The company is considering the weather condition to help predict the possibility of precipitations, which involves using various local climatological variables, including temperature, wind speed, humidity, dew point, and pressure. The data you will be handling was collected by a NOAA weather station located at the John F. Kennedy International Airport in Queens, New York.

Your task is to provide a high level analysis of weather data in JFK Airport. Your stakeholders want to understand the current and historical record of precipitations based on different variables. For now they are mainly interested in a macro-view of JFK Airport Weather, and how it relates to the possibility to rain because it will affect flight delays and etc.

Introduction

This project relates to the NOAA Weather Dataset - JFK Airport (New York). The original dataset contains 114,546 hourly observations of 12 local climatological variables (such as temperature and wind speed) collected at JFK airport. This dataset can be obtained for free from the IBM Developer Data Asset Exchange.

For this project, you will be using a subset dataset, which contains 5727 rows (about 5% or original rows) and 9 columns. The end goal will be to predict the precipitation using some of the available features. In this project, you will practice reading data files, preprocessing data, creating models, improving models and evaluating them to ultimately choose the best model.

Table of Contents:

Using this R notebook you will complete 10 tasks:

• 0. Import Modules
• 1. Download and Unzip NOAA Weather Dataset
• 2. Read Dataset into Project
• 3. Select Subset of Columns
• 4. Clean Up Columns
• 5. Convert Columns to Numerical Types
• 6. Rename Columns
• 7. Exploratory Data Analysis
• 8. Linear Regression
• 9. Improve the Model
• 10. Find Best Model

0. Import required modules

Tidymodels is a collection of packages that use tidyverse principles to easily do the entire modeling process from preprocessing initial data, to creating a model, to tunning hyperparameters. The tidymodels packages can be used to produce high quality statistical and machine learning models. Our Jupyter notebook platforms have a built-in Tidyverse, Tidymodels and rlang packages so we do not need to install these packages prior to loading library. However, if you decide to run this lab on your RStudio Desktop locally on your machine, you can remove the commented lines of code to install these packages before loading.

# Install tidymodels if you haven't done so
# install.packages("rlang")
# install.packages("tidymodels")

Note: After installing the packages, restart the kernel. Without installing the packages again, load them. Tidyverse and Tidymodels will be the two main packages you will use.

# Library for modeling
library(tidymodels)

# Load tidyverse
library(tidyverse)

── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──

✔ broom        1.0.7     ✔ recipes      1.1.0
✔ dials        1.3.0     ✔ rsample      1.2.1
✔ dplyr        1.1.4     ✔ tibble       3.2.1
✔ ggplot2      3.5.1     ✔ tidyr        1.3.1
✔ infer        1.0.7     ✔ tune         1.2.1
✔ modeldata    1.4.0     ✔ workflows    1.1.4
✔ parsnip      1.2.1     ✔ workflowsets 1.1.0
✔ purrr        1.0.2     ✔ yardstick    1.3.1

── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ purrr::discard() masks scales::discard()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ recipes::step()  masks stats::step()
• Use suppressPackageStartupMessages() to eliminate package startup messages

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.5
✔ lubridate 1.9.4     ✔ stringr   1.5.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ readr::col_factor() masks scales::col_factor()
✖ purrr::discard()    masks scales::discard()
✖ dplyr::filter()     masks stats::filter()
✖ stringr::fixed()    masks recipes::fixed()
✖ dplyr::lag()        masks stats::lag()
✖ readr::spec()       masks yardstick::spec()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Understand the Dataset

The original NOAA JFK dataset contains 114,546 hourly observations of various local climatological variables (including temperature, wind speed, humidity, dew point, and pressure).

In this project you will use a sample dataset, which is around 293 KB. Link to the sample dataset.

The sample contains 5727 rows (about 5% or original rows) and 9 columns, which are:

• DATE
• HOURLYDewPointTempF
• HOURLYRelativeHumidity
• HOURLYDRYBULBTEMPF
• HOURLYWETBULBTEMPF
• HOURLYPrecip
• HOURLYWindSpeed
• HOURLYSeaLevelPressure
• HOURLYStationPressure

The original dataset is much bigger. Feel free to explore the original dataset. Link to the original dataset.

For more information about the dataset, checkout the preview of NOAA Weather - JFK Airport.

1. Download NOAA Weather Dataset

Use the download.file() function to download the sample dataset from the URL below.

URL = 'https://dax-cdn.cdn.appdomain.cloud/dax-noaa-weather-data-jfk-airport/1.1.4/noaa-weather-sample-data.tar.gz'

url <- 'https://dax-cdn.cdn.appdomain.cloud/dax-noaa-weather-data-jfk-airport/1.1.4/noaa-weather-sample-data.tar.gz'
download.file(url, destfile = "noaa-weather-sample-data.tar.gz")

Untar the zipped file.

untar("noaa-weather-sample-data.tar.gz", tar = "internal")

Warning message in untar2(tarfile, files, list, exdir, restore_times):
“using pax extended headers”

2. Extract and Read into Project

We start by reading in the raw dataset. You should specify the file name as "noaa-weather-sample-data/jfk_weather_sample.csv".

noaa_data <- read_csv("noaa-weather-sample-data/jfk_weather_sample.csv")

Rows: 5727 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): HOURLYPrecip
dbl  (7): HOURLYDewPointTempF, HOURLYRelativeHumidity, HOURLYDRYBULBTEMPF, H...
dttm (1): DATE

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Next, display the first few rows of the dataframe.

head(noaa_data)

A tibble: 6 × 9

DATE	HOURLYDewPointTempF	HOURLYRelativeHumidity	HOURLYDRYBULBTEMPF	HOURLYWETBULBTEMPF	HOURLYPrecip	HOURLYWindSpeed	HOURLYSeaLevelPressure	HOURLYStationPressure
<dttm>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>	<dbl>	<dbl>	<dbl>
2015-07-25 13:51:00	60	46	83	68	0.00	13	30.01	29.99
2016-11-18 23:51:00	34	48	53	44	0.00	6	30.05	30.03
2013-01-06 08:51:00	33	89	36	35	0.00	13	30.14	30.12
2011-01-27 16:51:00	18	48	36	30	0.00	14	29.82	29.80
2015-01-03 12:16:00	27	61	39	34	T	11	NA	30.50
2013-02-15 20:51:00	35	79	41	38	0.00	6	29.94	29.92

Also, take a glimpse of the dataset to see the different column data types and make sure it is the correct subset dataset with about 5700 rows and 9 columns.

glimpse(noaa_data)

Rows: 5,727
Columns: 9
$ DATE                   <dttm> 2015-07-25 13:51:00, 2016-11-18 23:51:00, 2013…
$ HOURLYDewPointTempF    <dbl> 60, 34, 33, 18, 27, 35, 4, 14, 51, 71, 76, 19, …
$ HOURLYRelativeHumidity <dbl> 46, 48, 89, 48, 61, 79, 51, 65, 90, 94, 79, 37,…
$ HOURLYDRYBULBTEMPF     <dbl> 83, 53, 36, 36, 39, 41, 19, 24, 54, 73, 83, 44,…
$ HOURLYWETBULBTEMPF     <dbl> 68, 44, 35, 30, 34, 38, 15, 21, 52, 72, 78, 35,…
$ HOURLYPrecip           <chr> "0.00", "0.00", "0.00", "0.00", "T", "0.00", "0…
$ HOURLYWindSpeed        <dbl> 13, 6, 13, 14, 11, 6, 0, 11, 11, 5, 21, 7, 17, …
$ HOURLYSeaLevelPressure <dbl> 30.01, 30.05, 30.14, 29.82, NA, 29.94, 30.42, 3…
$ HOURLYStationPressure  <dbl> 29.99, 30.03, 30.12, 29.80, 30.50, 29.92, 30.40…

3. Select Subset of Columns

The end goal of this project will be to predict HOURLYprecip (precipitation) using a few other variables. Before you can do this, you first need to preprocess the dataset. Section 3 to section 6 focuses on preprocessing.

The first step in preprocessing is to select a subset of data columns and inspect the column types.

The key columns that we will explore in this project are:

• HOURLYRelativeHumidity
• HOURLYDRYBULBTEMPF
• HOURLYPrecip
• HOURLYWindSpeed
• HOURLYStationPressure

Data Glossary:

• 'HOURLYRelativeHumidity' is the relative humidity given to the nearest whole percentage.
• 'HOURLYDRYBULBTEMPF' is the dry-bulb temperature and is commonly used as the standard air temperature reported. It is given here in whole degrees Fahrenheit.
• 'HOURLYPrecip' is the amount of precipitation in inches to hundredths over the past hour. For certain automated stations, precipitation will be reported at sub-hourly intervals (e.g. every 15 or 20 minutes) as an accumulated amount of all precipitation within the preceding hour. A “T” indicates a trace amount of precipitation.
• 'HOURLYWindSpeed' is the speed of the wind at the time of observation given in miles per hour (mph).
• 'HOURLYStationPressure' is the atmospheric pressure observed at the station during the time of observation. Given in inches of Mercury (in Hg).

Select those five columns and store the modified dataframe as a new variable.

hrly_data <- noaa_data %>% 
                  select(HOURLYRelativeHumidity, HOURLYDRYBULBTEMPF, HOURLYPrecip, 
                         HOURLYWindSpeed, HOURLYStationPressure)

Show the first 10 rows of this new dataframe.

head(hrly_data, 10)

A tibble: 10 × 5

HOURLYRelativeHumidity	HOURLYDRYBULBTEMPF	HOURLYPrecip	HOURLYWindSpeed	HOURLYStationPressure
<dbl>	<dbl>	<chr>	<dbl>	<dbl>
46	83	0.00	13	29.99
48	53	0.00	6	30.03
89	36	0.00	13	30.12
48	36	0.00	14	29.80
61	39	T	11	30.50
79	41	0.00	6	29.92
51	19	0.00	0	30.40
65	24	0.00	11	30.35
90	54	0.06	11	30.03
94	73	NA	5	29.91

4. Clean Up Columns

From the dataframe preview above, we can see that the column HOURLYPrecip - which is the hourly measure of precipitation levels - contains both NA and T values. T specifies trace amounts of precipitation (meaning essentially no precipitation), while NA means not available, and is used to denote missing values. Additionally, some values also have "s" at the end of them, indicating that the precipitation was snow.

Inspect the unique values present in the column HOURLYPrecip (with unique(dataframe$column)) to see these values.

unique(hrly_data$HOURLYPrecip)

1. '0.00'
2. 'T'
3. '0.06'
4. NA
5. '0.03'
6. '0.02'
7. '0.08'
8. '0.01'
9. '0.07'
10. '0.16'
11. '0.09'
12. '0.22'
13. '0.02s'
14. '0.24'
15. '0.18'
16. '0.05'
17. '0.04'
18. '0.09s'
19. '0.11'
20. '0.14'
21. '0.25'
22. '0.10'
23. '0.01s'
24. '0.58'
25. '0.12'
26. '0.13'
27. '0.46'
28. '1.07'
29. '1.19'
30. '0.34'
31. '0.20'
32. '0.36s'
33. '0.42'
34. '0.17'
35. '0.27'
36. '0.35'
37. '0.31'
38. '0.33'
39. '0.23'
40. '0.26'
41. '0.28'
42. '0.75'
43. '0.19'
44. '0.36'
45. '0.03s'
46. '0.07s'
47. '0.54'
48. '0.59'
49. '0.21'

Having characters in values (like the "T" and "s" that you see in the unique values) will cause problems when you create a model because values for precipitation should be numerical. So you need to fix these values that have characters.

Now, for the column HOURLYPrecip:

1. Replace all the T values with "0.0" and
2. Remove "s" from values like "0.02s". In R, you can use the method str_remove(column, pattern = "s$") to remove the character "s" from the end of values. The "$" tells R to match to the end of values. The pattern is a regex pattern. Look at here for more information about regex and matching to strings in R.

Remember that you can use tidyverse's mutate() to update columns.

You can check your work by checking if unique values of HOURLYPrecip still contain any T or s. Store the modified dataframe as a new variable.

hr_data <- hrly_data %>%
  mutate(HOURLYPrecip = str_replace_all(HOURLYPrecip, "T", "0.0"),
         HOURLYPrecip = str_remove(HOURLYPrecip, "s$"))  
unique(hr_data$HOURLYPrecip)

1. '0.00'
2. '0.0'
3. '0.06'
4. NA
5. '0.03'
6. '0.02'
7. '0.08'
8. '0.01'
9. '0.07'
10. '0.16'
11. '0.09'
12. '0.22'
13. '0.24'
14. '0.18'
15. '0.05'
16. '0.04'
17. '0.11'
18. '0.14'
19. '0.25'
20. '0.10'
21. '0.58'
22. '0.12'
23. '0.13'
24. '0.46'
25. '1.07'
26. '1.19'
27. '0.34'
28. '0.20'
29. '0.36'
30. '0.42'
31. '0.17'
32. '0.27'
33. '0.35'
34. '0.31'
35. '0.33'
36. '0.23'
37. '0.26'
38. '0.28'
39. '0.75'
40. '0.19'
41. '0.54'
42. '0.59'
43. '0.21'

5. Convert Columns to Numerical Types

Now that you have removed the characters in the HOURLYPrecip column, you can safely covert the column to a numeric type.

First, check the types of the columns. You will notice that all are dbl (double or numeric) except for HOURLYPrecip, which is chr (character or string). Use the glimpse function from Tidyverse.

sapply(hr_data, typeof)

HOURLYRelativeHumidity

'double'

HOURLYDRYBULBTEMPF

'double'

HOURLYPrecip

'character'

HOURLYWindSpeed

'double'

HOURLYStationPressure

'double'

Convert HOURLYPrecip to the numeric type and store the cleaned dataframe as a new variable.

precip_data <- hr_data %>%
    mutate(HOURLYPrecip = as.double(HOURLYPrecip))

We can now see that all fields have numerical data type.

glimpse(precip_data)

Rows: 5,727
Columns: 5
$ HOURLYRelativeHumidity <dbl> 46, 48, 89, 48, 61, 79, 51, 65, 90, 94, 79, 37,…
$ HOURLYDRYBULBTEMPF     <dbl> 83, 53, 36, 36, 39, 41, 19, 24, 54, 73, 83, 44,…
$ HOURLYPrecip           <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,…
$ HOURLYWindSpeed        <dbl> 13, 6, 13, 14, 11, 6, 0, 11, 11, 5, 21, 7, 17, …
$ HOURLYStationPressure  <dbl> 29.99, 30.03, 30.12, 29.80, 30.50, 29.92, 30.40…

6. Rename Columns

Let's rename the following columns as:

• 'HOURLYRelativeHumidity' to 'relative_humidity'
• 'HOURLYDRYBULBTEMPF' to 'dry_bulb_temp_f'
• 'HOURLYPrecip' to 'precip'
• 'HOURLYWindSpeed' to 'wind_speed'
• 'HOURLYStationPressure' to 'station_pressure'

You can use dplyr::rename(). Then, store the final dataframe as a new variable.

names(precip_data) <- c('relative_humidity', 'dry_bulb_temp_f', 'precip', 'wind_speed', 'station_pressure')
head(precip_data)

A tibble: 6 × 5

relative_humidity	dry_bulb_temp_f	precip	wind_speed	station_pressure
<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
46	83	0	13	29.99
48	53	0	6	30.03
89	36	0	13	30.12
48	36	0	14	29.80
61	39	0	11	30.50
79	41	0	6	29.92

7. Exploratory Data Analysis

Now that you have finished preprocessing the dataset, you can can start exploring the columns more.

First, split the data into a training and testing set. Splitting a dataset is done randomly, so to have reproducible results set the seed = 1234. Also, use 80% of the data for training.

set.seed(1234)
precip_split <- initial_split(precip_data, prop = 0.8)
train_data <- training(precip_split)
test_data <- testing(precip_split)

Next, looking at just the training set, plot histograms or box plots of the variables (relative_humidity, dry_bulb_temp_f, precip, wind_speed, station_pressure) for an intial look of their distributions using tidyverse's ggplot. Leave the testing set as is because it is good practice to not see the testing set until evaluating the final model.

ggplot(gather(train_data, cols, value), aes(x = value)) + 
  geom_histogram(binwidth = 10, fill = "skyblue") +  # Change color here
  facet_grid(.~cols)

Warning message:
“Removed 1817 rows containing non-finite outside the scale range (`stat_bin()`).”

8. Linear Regression

After exploring the dataset more, you are now ready to start creating models to predict the precipitation (precip).

Create simple linear regression models where precip is the response variable and each of relative_humidity, dry_bulb_temp_f,wind_speed or station_pressure will be a predictor variable, e.g. precip ~ relative_humidity, precip ~ dry_bulb_temp_f, etc. for a total of four simple models. Additionally, visualize each simple model with a scatter plot.

lm_spec <- linear_reg() %>%
set_engine(engine = "lm")
train_fit <- lm_spec %>%
    fit(precip ~ relative_humidity, data = train_data)
train_fit
ggplot(train_data, aes(x=relative_humidity, y=precip)) +
    geom_point() +
    stat_smooth(method = "lm", col = "purple", se = FALSE)

parsnip model object


Call:
stats::lm(formula = precip ~ relative_humidity, data = data)

Coefficients:
      (Intercept)  relative_humidity  
       -0.0203017          0.0004074  

`geom_smooth()` using formula = 'y ~ x'
Warning message:
“Removed 1349 rows containing non-finite outside the scale range
(`stat_smooth()`).”
Warning message:
“Removed 1349 rows containing missing values or values outside the scale range
(`geom_point()`).”

lm_spec <- linear_reg() %>%
set_engine(engine = "lm")
train_fit <- lm_spec %>%
    fit(precip ~ dry_bulb_temp_f, data = train_data)
train_fit
ggplot(train_data, aes(x=dry_bulb_temp_f, y=precip)) +
    geom_point() +
    stat_smooth(method = "lm", col = "yellow", se = FALSE)

parsnip model object


Call:
stats::lm(formula = precip ~ dry_bulb_temp_f, data = data)

Coefficients:
    (Intercept)  dry_bulb_temp_f  
      4.752e-03        3.226e-05  

`geom_smooth()` using formula = 'y ~ x'
Warning message:
“Removed 1349 rows containing non-finite outside the scale range
(`stat_smooth()`).”
Warning message:
“Removed 1349 rows containing missing values or values outside the scale range
(`geom_point()`).”

lm_spec <- linear_reg() %>%
set_engine(engine = "lm")
train_fit <- lm_spec %>%
    fit(precip ~ wind_speed, data = train_data)
train_fit
ggplot(train_data, aes(x=wind_speed, y=precip)) +
    geom_point() +
    stat_smooth(method = "lm", col = "red", se = FALSE)

parsnip model object


Call:
stats::lm(formula = precip ~ wind_speed, data = data)

Coefficients:
(Intercept)   wind_speed  
  0.0022478    0.0003811  

`geom_smooth()` using formula = 'y ~ x'
Warning message:
“Removed 1349 rows containing non-finite outside the scale range
(`stat_smooth()`).”
Warning message:
“Removed 1349 rows containing missing values or values outside the scale range
(`geom_point()`).”

lm_spec <- linear_reg() %>%
set_engine(engine = "lm")
train_fit <- lm_spec %>%
    fit(precip ~ station_pressure, data = train_data)
train_fit
ggplot(train_data, aes(x=station_pressure, y=precip)) +
    geom_point() +
    stat_smooth(method = "lm", col = "blue", se = FALSE)

parsnip model object


Call:
stats::lm(formula = precip ~ station_pressure, data = data)

Coefficients:
     (Intercept)  station_pressure  
          0.6844           -0.0226  

`geom_smooth()` using formula = 'y ~ x'
Warning message:
“Removed 1353 rows containing non-finite outside the scale range
(`stat_smooth()`).”
Warning message:
“Removed 1353 rows containing missing values or values outside the scale range
(`geom_point()`).”

9. Improve the Model

Now, try improving the simple models you created in the previous section.

Create at least two more models, each model should use at least one of the different techniques:

1. Add more features/predictors
2. Add regularization (L1, L2 or a mix)
3. Add a polynomial component

Also, for each of the models you create, check the model performance using the training set and a metric like MSE, RMSE, or R-squared.

Consider using tidymodels if you choose to add regularization and tune lambda.

#poly
poly_percip <- lm_spec %>% fit(precip ~ poly(relative_humidity, 2, raw = TRUE), data = train_data)
res_train <- poly_percip %>% 
            predict(new_data = train_data) %>%
            mutate(truth = train_data$precip)
train_rmse <- rmse(res_train, truth = truth, estimate = .pred)
train_rmse
#mlr
mlr_train <- lm_spec %>% fit(precip ~ relative_humidity + dry_bulb_temp_f + wind_speed, data = train_data)
mlr_res <- mlr_train %>% 
            predict(new_data = train_data) %>%
            mutate(truth = train_data$precip)
mlr_rmse <- rmse(mlr_res, truth = truth, estimate = .pred)
mlr_rmse

A tibble: 1 × 3

.metric	.estimator	.estimate
<chr>	<chr>	<dbl>
rmse	standard	0.04223233

A tibble: 1 × 3

.metric	.estimator	.estimate
<chr>	<chr>	<dbl>
rmse	standard	0.04239945

10. Find Best Model

Compare the regression metrics of each model from section 9 to find the best model overall. To do this,

1. Evaluate the models on the testing set using at least one metric (like MSE, RMSE or R-squared).
2. After calculating the metrics on the testing set for each model, print them out in as a table to easily compare. You can use something like:

model_names <- c("model_1", "model_2", "model_3")
train_error <- c("model_1_value", "model_2_value", "model_3_value")
test_error <- c("model_1_value", "model_2_value", "model_3_value")
comparison_df <- data.frame(model_names, train_error, test_error)

3. Finally, from the comparison table you create, conclude which model performed the best.

# Polynomial Regression Model (Train)
poly_percip <- lm_spec %>% 
  fit(precip ~ poly(relative_humidity, 2, raw = TRUE), data = train_data)

res_train_poly <- poly_percip %>% 
  predict(new_data = train_data) %>%
  mutate(truth = train_data$precip)

train_rmse_poly <- rmse(res_train_poly, truth = truth, estimate = .pred)

# MLR(Train)
mlr_train <- lm_spec %>% 
  fit(precip ~ relative_humidity + dry_bulb_temp_f + wind_speed, data = train_data)

mlr_res_train <- mlr_train %>% 
  predict(new_data = train_data) %>%
  mutate(truth = train_data$precip)

train_rmse_mlr <- rmse(mlr_res_train, truth = truth, estimate = .pred)

# Polynomial Regression Model (Test)
res_test_poly <- poly_percip %>% 
  predict(new_data = test_data) %>%
  mutate(truth = test_data$precip)

test_rmse_poly <- rmse(res_test_poly, truth = truth, estimate = .pred)

# Multiple Linear Regression Model (Test)
res_test_mlr <- mlr_train %>% 
  predict(new_data = test_data) %>%
  mutate(truth = test_data$precip)

test_rmse_mlr <- rmse(res_test_mlr, truth = truth, estimate = .pred)

# Create a formatted tibble for Train RMSE results
train_comparison_df <- tibble(
  .metric = "rmse",
  .estimator = "standard",
  .estimate = c(train_rmse_poly$.estimate, train_rmse_mlr$.estimate)
)

# Print Train RMSE results in the exact format as your image
print(train_comparison_df)

# A tibble: 2 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard      0.0422
2 rmse    standard      0.0424

Author(s)

Yiwen Li

Contributions

Tiffany Zhu

© IBM Corporation. All rights reserved.